// bpcpp.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <cstdlib>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <string>
#include <algorithm>
#include <vector>

#include <Windows.h>

#include "tnt.h"
#include "jama_lu.h"
#include "jama_eig.h"
#include "jama_qr.h"
#include "jama_svd.h"
#include "jama_cholesky.h"

using namespace std; 
using namespace TNT; 
using namespace JAMA;

int M = 1663;  // velikost pole
int N = M;
int random=2;		// ma-li pocitat s nahodnou matici 0 1 2
					//0 vlastni
					//1 random
					//2 file

__int64 elapsedTime;
double timeR,improve;
LARGE_INTEGER cntStart, cntEnd;
LARGE_INTEGER fr;

//string cesta = "c:\\jama\\";
string cesta = "../";
string jazyk = "_cpp";
string pripona = ".csv";

double ran(){return((double)rand()/((double)(RAND_MAX)+(double)(1)));}

bool isSym(Matrix< double > XX){
bool Sym=1;
int c=XX.num_cols();
int r=XX.num_rows();
for (int i=0; i < c && Sym; i++){
	for (int j=0; j < r && Sym; j++){
		Sym=(XX[i][j]==XX[j][i]);		
	}
}
return Sym;
}

string int2string(const int& number)
{
   ostringstream oss;
   oss << number;
   return oss.str();
}

void vypisHard(Matrix< double > XX)
{	printf("\nMatrix");
	int c=XX.num_cols();
	int r=XX.num_rows();
	for (int i=0; i < c; i++){
		for (int j=0; j < r; j++){
			printf("\t%3.2f", XX[i][j]);
			if(j==(c-1)){printf("\n");}	   
		}
	
	}
}
void vypisA(Array2D< double > XX)
{//	if (random==0){
		printf("\n Array2D");
		int dim1=XX.dim1();
		int dim2=XX.dim1();
		for (int i=0; i < dim1; i++){
			for (int j=0; j < dim2; j++){
				printf("\t%3.2f", XX[i][j]);
				if(j==(dim1-1)){printf("\n");}
			}
		}	
//	}
}

void vypisA1(Array1D< double > XX)
{
	printf("\n Array1D");
	int dim1=XX.dim1();	
	for (int i=0; i < dim1; i++){	
		printf("\t%3.2f", XX[i]);			
	}
}	

void vypisV(Vector<double> XX){
if (random==0){
	printf("Vektor=(");	
	for (int i=1; i < XX.size()+1; i++){
		printf("%3.2f ", XX(i));
	}
	printf(")\n");
	}
}

void fileWrite(string nazew,Array1D< double > XX){
	
	string velikost = int2string(XX.dim1());
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open())
	{
		int dim1=XX.dim1();		
		for (int i=0; i < dim1; i++)
		{  
			myfile << XX[i];
			//myfile << ";" ;							
			myfile << "\n";					   
		
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

void fileWrite(string nazew,Array1D< int > XX){
	
	string velikost = int2string(XX.dim1());
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open())
	{
		int dim1=XX.dim1();		
		for (int i=0; i < dim1; i++)
		{  
			myfile << XX[i];
			//myfile << ";" ;							
			myfile << "\n";					   
		
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

void fileWrite(string nazew,Array2D< double > XX){
	
	string velikost = int2string(XX.dim1());
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open()){
		int c=XX.dim1();
		int r=XX.dim2();

		for (int i=0; i < c; i++){
			for (int j=0; j < r; j++){
				myfile << XX[i][j];
				myfile << ";" ;				
				if(j==(c-1))
				{
					myfile << "\n";
				}	   
			}
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}


void fileWriteDouble(string nazew,Array2D< double > XX){
	
	string velikost = int2string(XX.dim1()*2);
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open()){
		int c=XX.dim1();
		int r=XX.dim2();

		for (int i=0; i < c; i++){
			for (int j=0; j < r; j++){
				myfile << XX[i][j];
				myfile << ";" ;				
				if(j==(c-1))
				{
					myfile << "\n";
				}	   
			}
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}


void fileWrite(string nazew,Vector<double> XX){
	
	string velikost = int2string(XX.size());
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open())
	{
		int s=XX.size();		
		for (int i=0; i < s; i++)
		{
			myfile << XX[i];								
			myfile << "\n";	
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

void fileWrite(string nazew,Matrix <double> XX){

	string velikost = int2string(XX.num_rows());
	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open()){
		int c=XX.num_cols();
		int r=XX.num_rows();

		for (int i=0; i < c; i++){
			for (int j=0; j < r; j++){				
				myfile << XX[i][j]; //cerr <<  XX[i][j];				
				myfile << ";" ;				
				if(j==(c-1))
				{
					myfile << "\n";
				}	   
			}
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

void timeWrite(Vector<double> XX, int size){
	
	string nazew = "time";		
	string velikost = int2string(size);

	string fileN= cesta + nazew + jazyk + velikost + pripona;	
	char *nazev = (char *)fileN.c_str();
	
	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open())
	{
		int dim=XX.dim();		
		for (int i=1; i < dim; i++)
		{
			myfile << XX[i];								
			myfile << "\n";	
		}
		myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

void fileWrite(string nazew,double x,int size){
	
	string velikost = int2string(size);
	string fileN= cesta + nazew + jazyk + velikost + pripona;
	char *nazev = (char *)fileN.c_str();

	ofstream myfile (nazev);
	myfile.precision(17);
	if (myfile.is_open()){
			myfile << x;
			myfile.close();
	}
	else {
		cout << "Unable to open file";
		system("PAUSE");
	}
}

template <class T>
inline Array2D<T> toArray(const Matrix<T> &A)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();	
	Array2D<T> Arr(M,N);
	for (int i=0; i<M; i++)
		for (int j=0; j<N; j++)
			Arr[i][j] = A[i][j];
	return Arr;	
}

template <class T>
inline Matrix<T> toMatrix(const Array2D<T> &Arr)
{
    Subscript M = Arr.dim1();
    Subscript N = Arr.dim2();
	Matrix<T> R(M,N);
	//Array2D<T> Arr(M,N);
	for (int i=0; i<M; i++)
		for (int j=0; j<N; j++)
			R[i][j] = Arr[i][j];
	return R;
}

double round(double d, double pocetdc){
   const double m = pow(10, pocetdc);
   return (d < 0.0 ? ceil(d * m - 0.5) : floor(d * m + 0.5)) / m;
}

struct conversion_failure { };

template <typename T>
T from_string (const std::string & s)
{
	T result;
	std::istringstream stream (s);
	if (stream >> result) return result;
	throw conversion_failure ();
}

void Tokenize(const string& str,vector<string>& tokens,const string& delimiters = ";")
{
	string::size_type lastPos = str.find_first_not_of(delimiters, 0);	// Skip delimiters at beginning.
	string::size_type pos     = str.find_first_of(delimiters, lastPos);	// Find first "non-delimiter".

	while (string::npos != pos || string::npos != lastPos)
	{		
		tokens.push_back(str.substr(lastPos, pos - lastPos));			// Found a token, add it to the vector.
		lastPos = str.find_first_not_of(delimiters, pos);		
		pos = str.find_first_of(delimiters, lastPos);					// Find next "non-delimiter"
	}
}

Matrix< double> nacti (string nazev){
	// 	nacti config;	
	Matrix< double > XX;
	string line;
	int velikost;
		
	string fileN = cesta + "config.txt";
	char *nazevN = (char *)fileN.c_str();

	ifstream myfile (nazevN);	
	if (myfile.is_open())
	{	
		getline (myfile,line);
		getline (myfile,line);
		velikost = from_string <int> (line);
		
		myfile.close();
	}	
	else cout << "Unable to open file"; 
	XX.newsize(velikost,velikost);

	// nacteni matice o velikosti 

	ifstream indata;					
	
	string velikostN = int2string(velikost);
	fileN= cesta + "matice" +nazev+ velikostN + ".csv";
	char *nazew = (char *)fileN.c_str();

	cout << "Open file:" << nazew << endl;

	indata.open(nazew); // opens the file	
	if(!indata) { // file couldn't be opened
		cerr << "Error: file could not be opened" << endl;
		system("PAUSE");
		exit(1);
	}
	int i=0;
	indata >> line;
	//cerr.precision(20);
	while ( !indata.eof() ) { // keep reading until end-of-file				
		vector<string> tokens;
		Tokenize(line, tokens);
		for(int j = 0; j < (int) tokens.size(); j++){			
			XX[i][j] = from_string <double> (tokens[j]);			
		 }
		 indata >> line;
		 i++;
	}	
	indata.close();	
	return XX;	
}

void startT(){
	QueryPerformanceCounter(&cntStart); //zacatek mereni casu	
}
void stopT(){
	QueryPerformanceCounter(&cntEnd);  //konec mereni casu
	elapsedTime = (cntEnd.QuadPart - cntStart.QuadPart) * 1000000000 / fr.QuadPart;
	timeR = (double)elapsedTime/1000000. - improve ;	
}

int _tmain(int argc, _TCHAR* argv[]){
//%%%%%%%%%%%%priorita%%%%%%%%%%% Windows.h
SetPriorityClass(GetCurrentProcess(), HIGH_PRIORITY_CLASS);

BOOL result = QueryPerformanceFrequency(&fr);
if(result) {	
	printf("Hardware podporuje citac s vysokym rozlisenim. Frekvence citace: %lld\n",fr.QuadPart);	
  } else {
    printf("Hardware nepodporuje citac s vysokym rozlisenim.\n");
	system("PAUSE");
	return EXIT_SUCCESS;	
}  

//%%%%PRESNOST QueryPerformanceCounter
QueryPerformanceCounter(&cntStart); //zacatek mereni casu
QueryPerformanceCounter(&cntEnd);   //konec mereni casu
elapsedTime = (cntEnd.QuadPart - cntStart.QuadPart) * 1000000000 / fr.QuadPart;
improve = (double)elapsedTime/1000000. ;
startT();
stopT();
printf("Improve time: %3.6f ms \n", timeR);

int i;
int j;
Matrix< double > A(M,N);
Matrix< double > B(M,N);
Matrix< double > C;
Matrix< double > X;
Matrix< double > V;

Array1D<double> Re;
Array1D<double> Im;
Array2D<double> Ar;
Array2D<double> Br;
Array2D<double> Cr;
Array2D<double> Dr;

TNT::Vector<double> Vec;
TNT::Vector<double> Vb;
TNT::Vector<double> timeM;

timeM.newsize(27);
for (i=0; i < timeM.size(); i++)
	timeM[i]=-1;

printf("Start C++ benchmark\n");   

switch(random){
    case 0:		
		A.newsize(3,3);
		B.newsize(3,3);	
		
		A[0][0]=1;	A[0][1]=2;	A[0][2]=2;
		A[1][0]=4;	A[1][1]=5;	A[1][2]=6;
		A[2][0]=7;	A[2][1]=8;	A[2][2]=9;		

		B[0][0]=4;	B[0][1]=5;	B[0][2]=6;
		B[1][0]=2;	B[1][1]=8;	B[1][2]=9;
		B[2][0]=3;	B[2][1]=7;	B[2][2]=8;		
		
		M=3;
		N=M;
		C.newsize(M,N);
		X.newsize(M,N);
		break;		
		
    case 1:	
		startT();
		//A		
		for (i=0; i < M; i++)
			for (j=0; j < N; j++)
				A[i][j] = ran();
		stopT();		
		printf("Cas vytvoreni matice: %3.6f ms \n", timeR);		
   
	//B
		for (i=0; i < M; i++)
			for (j=0; j < N; j++)
				B[i][j] = ran(); 		
		C.newsize(M,N);

		break;
	case 2: 
		A=nacti("A");		//vypisHard(A);
		B=nacti("B");
		M=A.num_rows();
		N=M;

		C.newsize(M,N);
		X.newsize(M,N);
		
		break; 
}
//////start matic////////

//Generate
startT();
for (i=0; i < M; i++)
	for (j=0; j < N; j++)
		X[i][j] = ran();
stopT();
printf("Cas vytvoreni matice: %3.6f ms \n", timeR);
timeM[1]=timeR;

//1 na diag
int rows = X.num_rows();
int cols = X.num_cols();
startT();	
for (int i=0; i < rows; i++)
	for (j=0; j < cols; j++)
		(i==j ? X[i][j]=1 : X[i][j]=0);
stopT();
printf("Cas 1 na diagonale: %3.6f ms \n", timeR);
timeM[2]=timeR;
fileWrite("eye",X);//vypis(C);

//copy
startT();
	X=A;
stopT();
printf("Cas kopirovani matice: %3.6f ms \n", timeR);
timeM[3]=timeR;
fileWrite("copy",X);

startT();
	Ar=toArray(A);	
stopT();
printf("Cas 1 get array : %3.6f ms \n", timeR);
timeM[4]=timeR;
//fileWrite("array",Ar);

startT();
	Ar=toArray(A);
	Br=Ar.subarray(0,M/2-1,0,N/2-1);
	X=toMatrix(Ar);
stopT();
printf("Cas 1 get subarray : %3.6f ms \n", timeR);
timeM[5]=timeR;
fileWriteDouble("part",Br);	//vypis(X);

TNT::Matrix<double>* pZn;
pZn = new TNT::Matrix<double>(M,N);
for (i=0; i < M; i++)
	for (j=0; j < N; j++)
		(*pZn)[i][j] = ran();
startT();
	delete pZn;//Zn.newsize(0,0);//Zn.~Matrix();	//Zn.destroy();
stopT();	
printf("Cas zniceni matice: %3.6f ms \n", timeR);

startT();
	X=TNT::transpose(A);
stopT();
printf("Cas transponovani matice: %3.6f ms \n", timeR);
timeM[6]=timeR;
fileWrite("tra",X);	//vypis(X);

//norm1
startT();
	double f = 0;
	double s;
	for (int j = 0; j < N; j++){
		s = 0;
		for (int i = 0; i < M; i++){
			s += abs(A[i][j]);
		}
		f = max(f, s);
	}
stopT();
printf("Cas vlastniho vypoctu norm1: %3.6f ms \n", timeR);
timeM[7]=timeR;
fileWrite("norm1",f,A.num_cols());

//norm2
startT();
Ar=toArray(A);
JAMA::SVD<double> svdn(Ar);
	f=svdn.norm2();
stopT();
printf("Cas vlastniho vypoctu norm2: %3.6f ms \n", timeR);
timeM[8]=timeR;
fileWrite("norm2",f,A.num_cols());

//normInf
startT();
	f = 0;	
	for (int i = 0; i < M; i++){
		double s = 0;
		for (int j = 0; j < N; j++){
			s += abs(A[i][j]);
			}
		f = max(f, s);
	}
stopT();
printf("Cas vlastniho vypoctu normInf: %3.6f ms \n", timeR);
timeM[9]=timeR;
fileWrite("normi",f,A.num_cols());

//normF
startT();
	f=norm(A);
stopT();
printf("Cas vypoctu normFrob: %3.6f ms \n", timeR);
timeM[10]=timeR;
fileWrite("normf",f,A.num_cols());

//uminus
startT();
	for (int i = 0; i < M; i++)		
		for (int j = 0; j < N; j++)
			X[i][j] = -(A[i][j]);			
stopT();
printf("Cas unary minus vlastni: %3.6f ms \n", timeR);
timeM[11]=timeR;
fileWrite("uminus",X);//vypis(X);

//add
startT();
	X=TNT::add(A,B);
stopT();
printf("Cas scitani matic (add): %3.6f ms \n", timeR);
timeM[12]=timeR;
fileWrite("plus",X);

//minus
startT();
	X=TNT::minus(A,B);	
stopT();
printf("Cas odecitani matic: %3.6f ms \n", timeR);
timeM[13]=timeR; //vypis(X);
fileWrite("minus",X);

//X=A;
//TNT::mult_element_eq(X,B);//nasobeni ulozi do X
startT();
	X=TNT::mult_element(A,B);		
stopT();
printf("Cas nasobeni po prvku: %3.6f ms \n", timeR);
timeM[14]=timeR;
fileWrite("times",X);//vypis(X);


int atimes = 4;
startT();
	for (int i = 0; i < M; i++)		
		for (int j = 0; j < N; j++)
			X[i][j] = atimes*A[i][j];	
stopT();
printf("Cas nasobeni prvkem: %3.6f ms \n", timeR);
timeM[17]=timeR;//vypis(X);
fileWrite("multa",X);

//take mozno X=A*B;
startT();
	X=TNT::mult(A,B);
stopT();
printf("Cas nasobeni: %3.6f ms \n", timeR);
timeM[18]=timeR;//vypis(X);
fileWrite("mult",X);



startT();
Ar=toArray(A);
for (int i=0; i < rows; i++)
	for (j=0; j < cols; j++)
		(i==j ? X[i][j]=1 : X[i][j]=0);		
Br=toArray(X);

if(A.num_cols()==A.num_rows()){
	JAMA::LU<double> lus(Ar);
	Cr=lus.solve(Br);
}
else{
	JAMA::QR<double> qus(Ar);
	Cr=qus.solve(Br);
}
X=toMatrix(Cr);
stopT();
printf("Cas prevodu na inverzni matici: %3.6f ms \n", timeR);
fileWrite("inv",X);
timeM[19]=timeR;


Ar=toArray(A);
startT();
	JAMA::QR<double> qrr(Ar);
stopT();
if (qrr.isFullRank()){	
	fileWrite("rank",A.num_cols(),A.num_cols());
}
else{fileWrite("rank",-1,A.num_cols());
}
timeM[20]=timeR;


//DET
startT();
	Ar=toArray(A);
	JAMA::LU<double> ludet(Ar);	
		f = ludet.det();
stopT();
printf("Cas det: %3.6f ms \n", timeR);
timeM[21]=timeR;
fileWrite("det",f,A.num_cols());


//A(piv)=L*U  singular
printf("LU:\n");
X=nacti("LU");
Ar=toArray(X);
TNT::Array1D<int> pivot;
startT();
	JAMA::LU<double> lu(Ar);	
stopT();
if(lu.isNonsingular()==0){	
	Br = lu.getL();
	Cr = lu.getU();	
	pivot = lu.getPivot();
	X=toMatrix(Br);	
	V=toMatrix(Cr);
	printf("Cas LU: %3.6f ms \n", timeR);
	timeM[22]=timeR;	//vypis(X);	vypis(V);
	fileWrite("lu_l",X);
	fileWrite("lu_u",V);
	fileWrite("lu_p",pivot);
}
else{printf("LU nelze provest!\n");
}


//Q*R
printf("QR:\n");
Ar=toArray(A);
startT();
	JAMA::QR<double> qr(Ar);	
stopT();
Br = qr.getQ();
Cr = qr.getR();
X=toMatrix(Br);
V=toMatrix(Cr);
printf("Cas QR: %3.6f ms \n", timeR);
timeM[23]=timeR;
fileWrite("qr_q",X);//vypis(X);vypis(V);
fileWrite("qr_r",V);


X=nacti("CHOL");
printf("Cholesky:\n");
Ar=toArray(X);
startT();
	JAMA::Cholesky<double> chol(Ar);
stopT();
if (chol.is_spd()){
	Br = chol.getL();
	X=toMatrix(Br);	
	printf("Cas Cholesky: %3.6f ms \n", timeR);
	timeM[24]=timeR;	//vypis(X);
	fileWrite("chol",X);
}
else{printf("Cholesky nelze provest!\n");}

//SVD-alwayas exists U*S*V'
printf("SVD:\n");
Ar=toArray(A);
startT();
	JAMA::SVD<double> svd(Ar);
stopT();
svd.getS(Ar); //vrati matici		
fileWrite("svd_s",Ar);
svd.getU(Ar);
fileWrite("svd_u",Ar);
svd.getV(Ar); //vrati vektor
fileWrite("svd_v",Ar);
printf("Cas SVD: %3.6f ms \n", timeR);	
timeM[25]=timeR;	//vypis(X);

//if A=A' ...
//A=V*D*V    //V*V'=I
printf("EIG:\n");
//X=nacti("CHOL");
//if(isSym(A)==0){	//if(isSym(A)==0){ je spravne	
	Ar=toArray(A);	
	startT();
		JAMA::Eigenvalue<double> eig(Ar);	
	stopT();

	eig.getD(Ar);			
	fileWrite("eig_d",Ar);	
	eig.getV(Ar);	
	fileWrite("eig_v",Ar);

	eig.getImagEigenvalues(Im);
	eig.getRealEigenvalues(Re);
	fileWrite("eig_iv",Im);
	fileWrite("eig_rv",Re);

	printf("Cas Eigenvalue: %3.6f ms \n", timeR);	
	timeM[26]=timeR;		
//}else{printf("Eigenvalue nelze provest!\n");}
	
//ziskani vektoru
startT();
	Vec=A.diag();
stopT();
printf("Cas prevedeni diag na vektor: %3.6f ms \n", timeR);
fileWrite("diag",Vec);//vypisV(Vec);

timeWrite(timeM,A.num_cols());

printf("\n rows %i, cols %i \n",A.num_rows(),A.num_cols());
printf("============KONEC============\n");

system("PAUSE");

return EXIT_SUCCESS;
}
